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Abstract. In this paper, we present a multigrid l/-cycle preconditioner for the Hnear system 
arising from piecewise hnear nonconforming Crouzeix-Raviart discretization of second order 
elhptic problems with jump coefficients. The preconditioner uses standard conforming subspaces 
as coarse spaces. We showed that the convergence rate of the multigrid V-cycle algorithm wiU 
deteriorate rapidly due to large jumps in coefficient. However, the preconditioned system has 
only a fixed number of small eigenvalues, which are deteriorated due to the large jump in 
coefficient, and the effective condition number is bounded logarithmically with respect to the 
mesh size. As a result, the multigrid V-cycle preconditioned conjugate gradient algorithm 
converges nearly uniformly. Numerical tests show both robustness with respect to jumps in the 
coefficient and the mesh size. 



1. Introduction 

In this paper, we present a multigrid preconditioner for the linear system arising from piece- 
wise linear nonconforming Crouzeix-Raviart (CR) discretization of the following second order 
elliptic problems with jump coefficients: 

-V • (kVu) = fmn, u\9n = 0. (1.1) 

Here, Q C R*^ {d = 2,3) is an open polygonal domain and / G L'^{Q). We assume that the 
diffusion coefficient k G L°°(r2) is piecewise constant, namely, k{x)\q^ = Km is a constant for 
each (open) polygonal subdomain ^Im satisfying VJ^^-^Vlm = ^ and Q.m n = for m / n. 

Developing efficient multilevel solvers/preconditioners for the CR discretization is not only 
important for its own sake, but it has several important applications. In particular, it can be 



used to develop efficient solvers for other discretizations of (1.1). For example, by using the 



equivalence between CR discretization and mixed methods (cf. [Ij), the preconditioners for CR 



discretization can be applied to solve mixed finite element discretizations for (1.1) (see |lll I16j 
for the multigrid algorithms in the case of smooth coefficients) . This relationship is also used in 
|21| to analyze multigrid algorithms for mixed finite element discretization on adaptively refined 
meshes. Another important application is on the design of efficient multilevel solvers for interior 
penalty discontinuous Galerkin (IPDG) methods; we refer to [6j for the Poisson problem, and 



to 1^ for (1.1) with jump coefficients. 



A lot of work has been done to develop multilevel solvers and preconditioners for the piecewise 



linear CR discretization of (1.1) in the context of smooth (or constant) coefficients. Multigrid 
algorithms were developed and analyzed in [TOl 13 |9l [271 HI] ■, and two- level additive Schwarz 
preconditioners are presented in [l^. Hierarchical and BPX-type multilevel preconditioners 
were proposed in [22J, and their optimality was shown in [23j. Most of the aforementioned 
works define the multilevel structure using the natural sequences of nonconforming spaces. Since 
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these nonconforming spaces are non-nested, special intergrid transfer operators (restriction and 
profongation) are needed in the analysis and implementation of the algorithms. 

On the other hand, we observe that the standard piecewise linear conforming finite element 
space is a proper subspace of the CR space on the same triangulation. Therefore, one may be able 
to reduce the nonconforming discretization to the conforming one for which efficient and robust 
multilevel solvers/preconditioners are available. Due to the nested nature of the finite element 
spaces, one can also apply the standard multilevel analysis. Moreover, in the implementation of 
the algorithms one does not need special intergrid transfer operators; instead one could just use 
the natural inclusion operators. In this article, we propose a multigrid preconditioner for solving 
CR discretization, which consists of a standard smoother (Gauss-Seidel, or Jacobi iterative 
methods) on the fine nonconforming space and standard multigrid solver on the conforming 
subspaces. The only difference between this algorithm and the standard multigrid algorithm on 
conforming spaces is the prolongation and restriction operators on the finest level. Since the 
spaces are nested, the prolongation is simply the natural inclusion operator from the conforming 
space to the CR space. Some preliminary numerical results of this paper has been reported in [5] . 
We remark that the idea of using conforming subspace as preconditioner for CR discretization in 
the context of smooth coefficients has been used in [281 [301 [T71 [23] , aiid [25] for jump coefficient 
case. 

This article could be thought as an extension of the results in |3], where a BPX-type (additive) 
multilevel preconditioner for CR discretization was developed, and was applied to solve the whole 



family of interior penalty DG discretizations (both symmetric and nonsymmetric) for (1.1). No 
analysis on multiplicative multilevel preconditioner was given in [1]. In this article, we discuss 
the robustness of a multiplicative preconditioner, the multigrid V-cycle preconditioner. The 
analysis of the algorithm relies on some technical tools developed in [31j and [Jj. We show 
that the convergence rate of multigrid F-cycle algorithm will deteriorate rapidly due to large 
jump in the coefficient. On the other hand, we also show that the preconditioned system has 
only a fixed number of small eigenvalues deteriorated by the jump of the coefficients and mesh 
size, and the other eigenvalues are bounded nearly uniformly. Therefore, the resulting multigrid 
preconditioned PCG algorithm is robust with respect to the coefficient k and nearly optimal 
with respect to the mesh size h. 

Although the proof technique is similar to that in | 3Tj for the conforming case, we observe 



a major difference in the conclusion in Theorem 3.7 when d = 2. In the conforming case, the 
condition number of BA won't deteriorate due to the coefficients when d = 2, as was pointed 
out in [31, Remark 5.2]. However, the condition number of BA depends on the coefficient in 



the current nonconforming case, as we can see from the numerical tests (cf. Table 4.1 and 



Figure |4.2|). 

The remaining of this paper is organized as follows. In Section [2| we give basic notation, 
the finite element discretizations. We also review the PCG convergence theory, and some prop- 
erties of an interpolation operator introduced in [4j . In Section [sj we present the multigrid 
algorithm, discuss its implementation, and prove the convergence rate and the robustness as a 
preconditioner. Finally, we give some numerical experiments to verify our theory in Section [4] 

Throughout the paper, we will use the notation xi < yi, and X2 ^ y2, whenever there exist 
constants Ci, C2 independent of the mesh size h and the coefficient k or other parameters that 
^1, X2, yi and y2 may depend on, and such that xi < Ciyi and X2 > C2y2- 

2. Preliminaries 

In this section, we introduce some basic notation and finite element spaces. We give an 
overview of the PCG convergence theory. We also recall some approximation and stability 
properties of the interpolation operator introduced in [4J, which play an essential role in the 
convergence analysis of the multigrid algorithm. 
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Finite Element Spaces. The weak form of the model problem (1.1) reads: Find u G -ffg(J7) 
such that 

a{u,v) := {kVu,Vv) = {f,v) V^; € H^{n) . (2.1) 

We assume that there is an initial (quasi-uniform) triangulation To such that kt '■= h{x)\t is a 
constant for all T G To- Let 71 := 7/i, (Z = 1, • • • , L) be a family of uniform refinement of To with 
mesh size hi . Without loss of generality, we assume that the mesh size hi ~ 2^-^h (/ = O,-- - ,L) 
for h = hi- We denote £h the set of all edges (in 2D) or faces (in 3D) of Th- Let V^^ be the 
piecewise linear nonconforming Crouzeix-Raviart finite element space defined by: 

V^^=!^v E L'^in) : v^^ E pi(r) VT E % and jlvjeds = Ve E ^T^j , 

where P^(T) denotes the space of linear polynomials on T and {vje denotes the jump across 
the edge/face e ^ £h with |f]e = v when e C d^. Then the P^-nonconforming finite element 
discretization of (1.1) reads: 



Find u E ^ : ah{u, w) := / ktVu • Vw = (/, w), Vu; E ^. 



(2.2) 



The bilinear form a/i(-, •) induced a natural energy norm: := \/ ah{v, v) for any v E V^^. 

We denote the full (broken) weighted norm by ||f ||| ^ := ||?;||o ^ + |f || ^. For any w E H^{Q,), 
let ||t(^||o,K and Iil^Ii^k be the weighted norm and weighted semi-norm defined as 

Jn ' Jn 

Let A be the operator induced by the bilinear form ah{-,-), namely 

{Av, w) = {v, w)a ■■= ah{v, w), Vf , w E V^^. 

In the whole paper, we will use the notation || • m or | • 1^ interchangeably for the energy norm. 
Then solving (2.2) is equivalent to solve solve the linear system 

Au = f. (2.3) 

In the multigrid preconditioner defined in the next section, on each level / = 0, • • • , L we use 
standard conforming finite element discretization of (1.1) as coarse spaces: 

Find uiGVi: a{ui,vi) = if,vi), Vvi E Vi. (2.4) 

Here Vi := {v G Hq{Q) : v\t & Fi(T), VT E 71} is the standard continuous piecewise Pi- 
conforming finite element space defined on 71. For each / = 0, • • • , L, we define the induced 
operator for ( |2.4[ ) as 

{Aivi, wi) = a{vi,wi), yvi, wi E Vi. 

In the sequel, let us denote V^+i := for simplicity. We remark that all these finite element 
spaces are nested, that is, 

VoC-'-CVlCVl+i. (2.5) 

We also note that ah{vi, wi) = a{vi,wi) when vi,wi E V; for / = 0, 1, • • • , L since the space Vi is 
conforming. Therefore, with a little abuse of notation, we will use ah{-, •) as the bilinear form 
in all the finite element spaces including the conforming ones. 



4 Y. ZHU 

PCG Convergence. Let B be an symmetric positive definition (SPD) operator, known as a 



preconditioner. Instead of solving (2.3) directly, we consider solving the following preconditioned 



system by conjugate gradient method: 

BAu = Bf. 

Since A is SPD operator, BA is also SPD with respect to the energy inner product (•, ■)a- Based 
on the standard PCG theory, the convergence rate of PCG algorithm can be bounded as 

\\u-Uk\\A ^JlC{BA)-l \^ 
\\u-uo\\a~ \yjK:{BA) + 1 j ' 

where uq is the initial guess, {ufc} is the solution sequence, and K,{BA) is the (generalized) 
condition number of BA. 

This convergence rate is not sharp, in particular, when BA contains some extreme eigenvalues. 
Suppose that the spectrum a{BA) of BA, is divided in two sets: a{BA) = aQ{BA) U (Ji{BA), 
where aQ{BA) = {Ai, . . . , Am} contains of all "bad" eigenvalues, and ai{BA) = {Am+i, . . . , Xn} 
contains the remaining eigenvalues, which are bounded below and above uniformly. Namely, 
there exist some constants < a < 6 such that Xj S [a, h] for j = m + 1, . . . , A^, with N = 
dim{V^^). Then, the error at the k-th. iteration of PCG algorithm is bounded by (see e.g. [21 

/ I \ k—m 

\\u - u^Wa < 2{}C{BA) - ir 1 ~ ^ I lln-^xolU- (2.6) 



This convergence rate estimate ( |2.6| ) implies that if there are only a few small eigenvalues of 
BA in ao{BA), then the asymptotic convergence rate of the resulting PCG method will be 

dominated by the factor i,e., by \fhja where h = XjqiBA) and a = Xm^i{BA). We 

define (b/a), which determines the asymptotic convergence rate, as effective condition number 

(cf. m)-- 



Definition 2.1. Let y be a real A^-dimensional Hilbert space, and T : V ^ V he w. SPD 
linear operator, with eigenvalues < Ai < • • • < Aat. The m-th effective condition number of T 
is defined by 

Xn{T) 



)Cm{T) :-- 



Xm+l{T) 



Based on the estimate in (2.6 ), we need to analyze the spectrum of the preconditioned system 



carefully. That is what we are going to do for multigrid preconditioner. 

An Interpolation Operator. Given any / = 0, 1, • • • , L, let us review some main properties 
of an interpolation operator Qf : — >• Vi introduced in [3], which play a key role in the 
analysis of multilevel preconditioner s. 

Lemma 2.2. There exists an interpolation operator Qf : — )• Vi satisfying the following 
approximation and stability properties: 

Approximation: - Q^)v\\o,^, < hi\log{2hi/h)\^\v\\h,K, Vt; E Vf ^, (2.7) 

Stability: \Qfv\i,,<\\ogi2hi/h)\l\\v\\h,^ WveV^^, (2.8) 

with constants independent of the coefficient k and mesh size. 

A construction of such an operator and the proof of the above results are given in the Appendix 
of [1]. We would like to point out that the operators (/ = 0, 1, • • • , L) are not used in the 
actual implementation of the preconditioner. 
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Observe that on the right hand side of (2.7) and (|2.8|), the bounds are given in terms of the 



(broken) weighted full if^-norm In general, one cannot replace the norm by the 

energy norm |f I^^k induced by the bilinear form o/i(-, •). To replace this full weighted norm by 
the energy norm, we may use the Poincare-Priedrichs inequality for the nonconforming finite 
element space (cf. IIS [13]) to obtain: 



I l|2 ^ ( \ I \ \'^A < ( ^ IIT7 l|2 < maXTt^Thf^T, |2 

k; < maxKy / \v\ ax <- I max kt ] > Vt; L ^ < — : \v\h i^- 



By this inequality, we have 

Corollary 2.3. There exists an interpolation operator Qf : — )• Vi satisfying the following 
approximation and stability properties: 

\\{I-Qf)v\\o,. < jHK)hi\log{2hi/h)\'2\v\h,^ , yy(zv^R, (2.9) 

\Qfv\i,n < jH'^)nog{2hi/h)\"2\v\h,n , yveV^"^, (2.10) 

where = maxj^gT^ kt/ minTGTh '^T is the jump of the coefficient. 



yen 



Note the approximation and stability properties in Corollary 2.3 depend on the jump J^{k) 
of the coefficient. On the other hand, we may impose certain constraints on the finite element 
space to get robust approximation and stability properties. For this purpose, let I be the 
index set of floating subdomains (cf. [26] ) , namely the subdomains not touching the Dirichlet 
boundary: 

I:={i : measd-i{dnndni) = 0} . (2.11) 
We introduce the subspace V^'^ C V^^: 

\vG Vf ^ : / vdx = 0, yi£l\ , (2.12) 

The key feature of this subspace ( |2.12 ) is that the Poincare-Friedrichs inequality for the non- 
conforming finite element space (cf. |18| [T3] ) holds on each subdomain, and this allows us to 
replace the full weighted i?^-norm ||f by the energy norm |f for any v S V^^. 

We remark that the condition vdx = in (2.12) is not essential; other conditions could be 
used (see [SS]) as long as they allow for the application of a Poincare-type inequality. At this 
point, we would like to emphasize that the dimension of is related to the number of floating 
subdomains and in fact, dim(y^^) = dim(V^'^) — mo, where mo = f^T is the cardinality of X. 

By restricting now the action of the operator QJ^ to functions in V^^, we have the following 
result, as an easy corollary from Lemma |2.2[ 

Corollary 2.4. Let C be the subspace defined in (2.12). There exists an operator 

Ql '■ ~^ ^ satisfying the following approximation and stability properties: 

||(/-Qr)^^||o,. < hi\log{2hi/h)\'2\v\h,. , yveV^"^, (2.13) 
mv\i,^ < \log{2hi/h)\'2\v\h,. , yveV^"^. (2.14) 

Proof. By the deflnition of V^^, any v G ^h '^ satisfies the Poincare-Friedrichs inequality on 
each subdomain (cf. [T8l [T3]): 

/ \v\'^dx< I \Vv\^dx, i = l,---,M. 

Since k is piecewise constant on each subdomain, we have 

\\v\\o,^ < \v\h,^, yv e V^^. 
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The conclusion then follows by Lemma 2.2 and the above inequality. □ 



3. A MULTIGRID PrECONDITIONER 

The action of standard multigrid IZ-cycle iterative operator B := -Bl+i : Vl+i ^ Vl+i on a 
given g G Vl+i can be defined recursively by the following algorithm (cf. [8]): 

Algorithm 3.1 (F-cycle). Let gi+i = g, and Bq = . For / = 1, • • • , L + 1, we define 
recursively Bigi for any gi £ Vi by the following three steps: 

(i) Pre-smoothing : wi = Rigr, 

(ii) Subspace correction: W2 = wi + Bi-iQi-i{gi — Aiwi); 

(iii) Post-smoothing: Bigi := W2 + RJigi — AiW2)- 

Here, T* is the adjoint operator of T with respect to the energy inner product o/i(-, •), that 

is, 

ah{T*v,w) := ah{v,Tw), Mv,w £ V^^. 
In this algorithm, Qi : L^(J7) — ?• Vi is the standard projection defined by: 

{Qiv,wi) = {v,wi), \/wieVi, (/ = 0,---,L). 

The operator i?/ : — >• for / = 0,-- - ,L + lis called a smoother. Here we consider Ri to be the 
Gauss-Seidel smoother, which can be interpreted as an successive subspace correction algorithm 
based on a one-dimensional subspace decomposition = Y17=i with T/^* = span{0j}"i^^ (cf. 
[29]). Here n/ denotes the dimension of the space Vi. 



Algorithm 3.2 (Gauss-Seidel). Let wq = £ Vi. Then for any gi £ Vi, we define Rigi 
as the result of the following loop: for i = 1, • • • ,ni do 

(i) Find Cj G : a/,(ei,(/)J) = {gi,(t>l) - atiwi-i, 4>l); 

(ii) Update: Wi := Wi-i + e,. 

For i = 1, • • • , n;, let P/ : Vl+i V^ be the orthogonal projection defined by 

Using this projection operator, the solution Cj to the first step in Algorithm 3.2 is 

ei = Pi{Al^g-Wi.i). 



Hence, the following recursive relation holds from Algorithm 3.2 



A ^91 -Wi = {I - Pi){Ai ^gi - wi^i), i = 1, • • • ,nz. 



Therefore, the Gauss-Seidel smoother Ri -.Vi defined in Algorithm 3.2 can be written as 

Ri = [l - Pi)^ Ai\ (3.1) 

On each level Z = 0, 1, • • • , L -|- 1, we define the Galerkin projection Pi : Vl+i — ?• as 

ah{Piu,vi) = ah{u,vi), Mvi G Vi. 

The implementation of Algorithm |3.1| is almost identical to the implementation of standard 
multigrid T^-cycle (cf. [E]). We use standard prolongation and restriction matrices for con- 
forming finite element spaces Vi for I = 0, • • • , L. The corresponding matrices between Vl 
and Vl+1, are however different. The prolongation matrix on Vl can be viewed as the matrix 
representation of the natural inclusion : — )• Vl+i, which is defined by 

{Zlv){x) = ^ u(me)-0e(2;). 
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where ipe is the CR basis on the edge/face e & £h and nie is the barycener of e. Therefore, 
the profongation matrix has the same sparsity pattern as the edge-to- vertex (in 2D), or face-to- 
vertex (in 3D) connectivity, and each nonzero entry in this matrix equals the constant 1/d where 
d is the space dimension. The restriction matrix is simply the transpose of the prolongation 
matrix. 

Now we are in position to discuss convergence of the multigrid y-cycle algorithm. The 
analysis is based on the following identity [32j. 



subspaces satisfying W = Ylj=o ^j- Pj : W ^ Wj be the orthogonal projection with respect 



Lemma 3.3. Let W be a Hilbert space, and Wj C W for j = 0, ■ ■ ■ , J be a number of closed 

.,=o^i- LetPr-W 
to the inner product of W. Then 

\\{I-Pj)---{I 

where 



P 



0)\\C(W) 



1 



1 

1 + co' 



J 

Co = sup inf 

llf|l=i''=E/=o''j j=o 




Recall that from (2.5), the finite element spaces {V/}^^ are nested. Then we have the 
following space decomposition (recall that Vl+i 



CR 



L+1 
1=0 



(3.2) 



rithm 3.1 has the representation (cf. [HI (3.4)]) 



Combining with (3.1) for the smoother Ri, the multigrid l^-cycle operator B := i?L+i in Algo- 

(3.3) 



/ L+1 ni \ * / L+1 ni \ 

I - BA = ui - p,)w\{{i - pi)\ (/-Po)nn(^-^/) ' 

V 1=1 i=l J \ 1=1 1=1 J 

and the error propagation of the multigrid V-cycle iteration has the form: 

u-Uk = {I - BA)^{u - uo). 

We define the convergence rate p := ||/ — i^^m, then obviously the above iteration is convergent 
if /9 < 1. From (3.3), we have the following simple relationship between the extreme eigenvalues 
of BA with the convergence rate p: 

Proposition 3.4. We have the following estimate for the maximum and minimum eigenvalues 
ofBA: 

Amax(S^) < 1, and XmUBA) = l-p. (3.4) 
Proof. Obviously, for any v G we have 

2 



Ohiil- BA)v,v) 



' L+1 ni \ 

y 1=1 i=l / 



> 0. 



Therefore, ah{BAv,v) < ah{v,v), Vv G Vh^, which implies that Amax(-B^) < 1. 
By definition of ||/ — -B^||a, we have 



\I -BA\ 



A = sup 

O^veV^R 



ahijl - BA)v,v) 
ah{v,v) 



inf 



ahiBAv,v) 
ah{v,v) 
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From this identity, it is obvious that 

Amm(-B^) ■ 

This completes the proof. 



. ah{BAv,v) 

mf — = 1 - p. 

o^vGV^'^ ah{v,v) 



□ 



On the other hand, notice that I — BA is self-adjoint with respect to the inner product a/i(-, •), 
we have 

2 



\I-BA\\^= sup ah{{I - BA)v,v) 
\Ma=^ 



L+l ni 
1=1 i=l 



From Lemma 3.3, we deduce that 



p= \\I -BA\ 



Co 



^ — , with Co = sup inf c(t'), 



(3.5) 



where 



L+l ni 



c{v) = mv-vo)\i, + Y,Y. 



=1 i=i 



E 



,(fc,i)>a,i) 



h,K 



with € span{(/);[}. In other words, we only need to estimate the upper bound of cq to obtain 



an upper bound of convergence rate p of the multigrid F-cycle algorithm. By Proposition 3.4 
the upper bound of p also provides a lower bound of minimum eigenvalue of BA. 

For this purpose, we introduce the following lemma to bound c{v). To simplify the notation, 
we set Q2+1 ~ ^'^^ Q-i ~ 0- Then for any v G Vl+i, we can decompose it as 



L+l L+l 

^ = Y.^Qi -Qti)v = ^vu where vi = {Qf - Qti)v ■ 



(3.6) 



1=0 1=0 
Clearly, G for / = 1, • • • , (L + 1) and vq = QqV G Vq. 



Lemma 3.5. Given any v G Vl+i with the decomposition (3.6), we have 



L+l L+l 

iv) < \Piv - Q^v\i, + hf II (Qr - Qt-Mt^ ■ 

1=0 1=1 



(3.7) 



Proof. By the decomposition (3.6), we can write v G as 

L+l L+l ni 

w = J]] = WO + XI 

1=0 1=1 i=l 

with vi = {Qf — Qi_^)v. Noticing that 



L+l Uk 

E <--=T.Yl < = - + E 

(k,j)>{l,i) k=l j=i+l j=i+l 
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we have 



L+l ni 

EE 

1=1 i=l 



Pi E 



{k,j)>{i,i) 



h,K 



L+l ni 

EE 

1=1 1=1 



Pi E 



L+l ni 

S^EE 

1=1 i=l 



\Pii--Q»\K.n,+ 



ik,j)>ii,i) 



h,K,Ql ; 



Pi E 



j=i+l 



L+l ni 

<2EEI 



+ 



1=1 i=l 



E-/ 

j=i+i 



h,K,Q,i i . 



where ^i^i := supp(0j) and \wW ^ , := X^TeTl rcH; ^ /r '^t| Vwpdx. By quasi-uniformity of the 
triangulations, it is not difficult to verify that 



Y.\Pi^^-Qi<,.,n,.^\Pi^-Qi<^ 

1=1 

and by inverse inequahty and quasi-uniformity of the triangulations 

2 



E 

1=1 



ni 

E^ 

j=i+i 



h,K,fti . 



ni ni 

^ E E 

i=l j=i+l 



<r- i,-2 II ||2 



hf\\{Qt-Qti>l^^. 



Therefore, we get 



L+l 



L+l 



iv) < E 1^'^ - Q>\L + E V' UQf - Qi-M\o,. ■ 



1=0 



1=1 



This completes the proof. 



□ 



Based on this lemma, in order to estimate the convergence rate of the multigrid algorithm, 
we only need the stability and approximation properties of the interpolation operators Qf 
(/ = 0, • • • , L). Now we are in position to give an estimate on the convergence rate p of the 
multigrid F-cycle Algorithm |3.1[ 



Theorem 3.6. Let B he the multigrid V -cycle iterator defined in Algorithm \3.1\ Then the 
convergence rate p = \\I — BA\\a satisfies 



1 + CoJ{k)L^'' 



(3.8) 



where Co is a constant independent of the coefficient k and the mesh size h (or number of levels 
L). Moreover, the condition number IC{B A) satisfies 

]C{BA) < 1 + CoJ{k)L^. (3.9) 



Proof. Given any v G V^^, we make use of Corollary 2.3 to bound the right hand side of ( |3.7 
For the first term, by the stability (2.10) of Qf we have, for / = 0, 1, • • • ,L 



\Piv - Qfvl < \vl + \Qfv\, < 1 + J5 (k) |log hi\-^ \vl 
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Notice that Q^+i ~ obtain that 

L+l L 

J2\Piv-Qfv\l, = Y,\Piv-Qfv\l. + \PL+i 



V — V 



h,K 



l=Q 



1=0 



< J{^)(^J2\^oghi\^\v\l,<J{^)L' 



\h,K 



For the second term in the right hand side of (3.7), by the approximation property (2.9) of 

\\Qfv - Ql,v\\^^^^ < \\v-Qtv\\f^^, + \\v-Qt_^v\l^^ 

< j5(K)/lj_i|log/l,_i| 



1 

|2 \v\ 



h,K ■ 



Hence we have 



L+l /L+l \ 

i=\ \i=\ ) 



r2 I |2 



Therefore, we obtain cq < CqJ{k)1? for some constant Co independent of coefficient k and 
meshsize h. The estimate ( |3.8| ) then follows by the identity (3.5). Finally, the condition number 
estimate follows by (3.4) in Proposition 3.4 and (3.8). □ 



From Theorem |3.6[ we conclude that the convergence rate of the multigrid V-cycle iteration 
is deteriorated by the jump of the coefficients rapidly. This phenomenon is justified also by the 
numerical experiment in Section [4j On the other hand, we can show that this deterioration due 
to the large jump in coefficient is only limited to a few fixed number of small eigenvalues in 
the preconditioned system. Therefore, the asymptotic convergence rate of the PCG algorithm 
is nearly uniform, as stated in the following theorem. 



Theorem 3.7. Let B be the multilevel preconditioner defined by Algorithm \3.1\ The effective 
condition number K-mo {BA) satisfies 

lCm,{BA)<ClL^, (3.10) 

where mo = #2^ is the number of floating subdomains (cf. ( 2.11[ ) ). Therefore, the PCG algorithm 
has the following convergence rate estimate: 

II II / r 1 \ k—mo 

"^-^^"^<2(CoJ(^)L2r«^^^^"^ 



\u - Uo\\A 



CiL + l 



Proof. To obtain the effective condition number, we restrict ourself in the subspace C 
defined by (2.12). By the identity (3.5), we have 

2 



\I -BA\ 



L+l ni 

{i-p')\{\{{i-pi 

1=1 i=l 



Co 



l + co' 



with 



Co 



sup inf 



c(v 



Given any v G V^^, we make use of Corollary 2.4 to bound the right hand side of (3.7). For 
the first term, by the stability (2.14) of Qf we have, for / = 0, 1, • 



L 



\Piv - Qfv\^,, < \v\^,^ + igr^i,,. < 1 + \ioghi\-^\ \v\^^^ 
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Notice that Ql^i = I, we obtain that 

L+l L 

J2\P,v-Qfv\l^ = J2\Plv-Qfv\l. + \PL+i 



V — V 



h,K 



1=0 



< 



j=o 




For the second term in the right hand side of (3.7), by the approximation property (2.13) of 
Ql, (/ = !,••• + 

||Q^-Qr-i^llo,. < \\^-Qf^\\o,K + h-Qt-Ao,. 



< hi_i |log/ii_i|2 \v\ 



Hence we have 



L+l /L+l \ 

hf II (Qf - QtiHl. ^ E I log^'-ii 
1=1 \i=i J 



|2 ^ r2 I |2 



h,K ' 



Therefore, we obtain that cq < L^. Then by mini-max Theorem (cf. [191 Theorem 8.1.2]) we 
obtain 

/ ^ a(BAv, v) 1 ^ ^ 9 

o^v&vc^R a{v, V) 1 + Co 



Together with (3.4), we get the estimate on the effective condition number of BA: 

lC^,{BA)<ClL\ 

for some constant Ci independent of the coefficient k and mesh size h (or levels L). The con- 
vergence rate estimate of the PCG algorithm follows by (3.9), (3.10) and ( |2.6[ ). This completes 
the proof. □ 

Remark 3.8. We remark that mo = #X gives an upper bound of the total number small eigen- 
value. In some cases, the number of bad eigenvalues might be less than mo, as one can observe 
from the numerical tests. Since mo is a fixed integer depending only on the distribution of 
the coefficient, the asymptotic convergence rate of the PCG algorithm is dominated by , 
which is independent of the coefficient. 

4. Numerical Results 
In this section, we present some numerical experiments in both 2D and 3D domains to justify 



the performance of the multigrid F-cycle preconditioner in Algorithm 3.1 



In the first example, we consider the model problem (1.1) in the square 
coefficients k = 1 in the subdomain [—0.5, 1]^ U [0, 0.5]^, and n 



[-1,1]2 with 
£ for the remaining subdomain, 



where e varies (see Figure 4.1). We start with a structured initial triangulation on level 
with mesh size h = which resolves the jump interface of the coefficient. Then on each 
level, we uniformly refine the mesh by subdividing each element into four congruent ones. In 
this example, we use 1 forward/backward Gauss-Seidel iteration as pre/post smoother in the 
multigrid preconditioner, and the stopping criteria of the PCG algorithm is ||?"A:||/lko|| < 10^^ 
where = f — Auk is the the residual at A;-th iteration. 



Figure 4.2 shows the eigenvalue distribution of the multigrid y-cycle preconditioned system 
BA when h = and e = 10~^. As we can see from this figure, there is only one small eigenvalue 
deteriorated by the coefficient and mesh size. This is different from the conforming case. As 
we know in the conforming case (cf. [31', Remark 5.2]), there will be no such extremely small 
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eigenvalue. However, we still observe a small eigenvalue in Figure 4.2 in this nonconforming 
case. 




1 1 

0.9 - 
0.8 
0.7 - 
0.6 - 

0.5 * ■ 
0.4 - 
0.3 
0.2 - 
0.1 
0M« 



10 



15 



20 



Figure 4.1. 2D Computa- 
tional Domain 



Figure 4.2. Eigenvalue 
Distribution of BA in 2D 



Table 4.1 shows the estimated condition number /C (number of PCG iterations) and the 
effective condition number /Ci of the multigrid preconditioned system on V^^. As we can see 
from this table, although the condition number is deteriorated due to the jumps in the coefficient 
and mesh size, the corresponding effective condition number is nearly uniform, which coincides 
with the theory. 



e 


levels 





1 


2 


3 


4 


1 


/C 

/Ci 


1.65 (8) 
1.44 


1.83 (10) 
1.78 


1.9 (10) 
1.77 


1.9 (10) 
1.78 


1.89 (10) 
1.76 


10-1 


/C 
/Ci 


3.78 (10) 
1.89 


3.69 (11) 
1.87 


3.76 (12) 
1.93 


3.79 (12) 
1.92 


3.88 (12) 
1.95 


10-2 


/C 
/Ci 


23.4 (12) 
2.15 


23.6 (13) 
1.96 


24.6 (13) 
1.99 


25.1 (14) 
1.97 


26 (15) 
2.24 


10-3 


/c 


218 (13) 
2.19 


223 (14) 
1.98 


232 (15) 
2 


238 (16) 
1.98 


246 (16) 
2.29 


10-4 


/c 


2.17e+03 (14) 
2.2 


2.21e+03 (15) 
1.98 


2.31e+03 (16) 
2 


2.37e+03 (18) 
1.98 


2.45e+03 (18) 
2.3 


10-5 


/c 


2.17e+04 (15) 
2.2 


2.21e+04 (16) 
1.98 


2.31e+04 (17) 
2 


2.37e+04 (20) 
1.98 


2.76e+04 (21) 
2.64 



Table 4.1. 
fCi in 2D 



Estimated condition number /C and the effective condition number 



As a second example, we consider the model problem in a 3D unit cube with similar setting. 
We set the coefficient k = 1.0 in subdomains Jli = [0.25, 0.5]^ and 0,2 = [0-5, 0.75]^; and k = e for 
the remaining subdomain (see Figure 4.3). The initial mesh (level =0) has a mesh size ho = 2~'^, 
which resolves the jump interface. In this example, we used 5 forward/backward Gauss-Seidel 
as smoother in the multigrid preconditioner, and the stopping criteria ||?^fc||/lko|| < IQ-^^ for 
the PCG algorithm. 
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Figure 4.3. 3D Computa- 
tional Domain Figure 4.4. Eigenvalue 

Distribution of BA in 3D 



e 


levels 





1 


2 


3 


1 


/C 
ICi 


1.19 (8) 
1.16 


1.34 (11) 
1.26 


1.37 (11) 
1.31 


1.36 (11) 
1.29 


10-1 


K 
fCi 


2.3 (10) 
1.60 


1.94(13) 
1.56 


1.75 (13) 
1.45 


1.67 (14) 
1.43 


10-^ 


IC 
fCi 


86.01 (11) 
2.4 


63.07 (16) 
2.12 


52.67 (17) 
1.89 


48.19(17) 
1.78 


10-5 


IC 
ICi 


8.39+03 (13) 
2.44 


6.15e+03 (18) 
2.14 


5.13e+03 (19) 
1.91 


4.70e+03(19) 
1.80 


10-7 


IC 
ICi 


8.39+05 (14) 
2.45 


6.15e+05 (21) 
2.14 


5.13e+05 (23) 
1.91 


4.70e+05(21) 
1.80 



Table 4.2. Estimated condition number IC (PCG iterations) and effective con- 
dition number ICi for multigrid y-cycle preconditioner in 3D 



Figure 4.4 shows the eigenvalue distribution of the multigrid T^-cycle preconditioned system 
BA when h = and e = 10"^. As before, this figure shows that there is only one small 
eigenvalue deteriorated by the coefficient and mesh size. 

Table 4.2 shows the estimated condition number IC (with the number of PCG iterations), and 
the effective condition number ICi. As we can see from this table, the condition number IC gets 
bigger as the jump J{k) = 1/e becomes larger. Thus the condition number deteriorates due to 
the jump of coefficient. However, the effective condition number fCi remains in a small range. 
This justifies that the PCG algorithm with multigrid y-cycle preconditioner is a robust solver 
for solving (1.1) with jump coefficients. 

Table 4.3 and 4.4 shows the estimated convergence rate of the multigrid l^-cycle with 2 Gauss- 
Seidel smoothers and 5 Gauss-Seidel smoothers respectively, with respect to different coefficients 
e = 10~*, i = 0, 1, • • • ,5 (rows) and different levels L = 0, 1, • • • ,3 (columns). The convergence 
rate of the multigrid F-cycle deteriorates rapidly with respect to the jump. Actually, when 
e < IQ-^ the multigrid V-cyc\e is unfavorable. When the jump is not so severe, the more 
smoothing step, the faster the algorithm converges. However, for problems with large jump, 
more smoothing steps does not improve the convergence rate much. 
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Another interesting phenomenon we observe from this two tables is that the convergence rate 
seems to be uniform with respect to the number of levels, when the jump is fixed. The same 
phenomenon can also be observed from Table |4.2| for the effective condition numbers. This is 
also different from the conforming case in [3H Table 2], where we could observe some growth of 
the convergence rate with respect to the number of levels. The theoretical investigation of this 
phenomenon will be a future work. 







e 


Levels 


h 


1 


10-^ 


10-^ 


10-3 


10-4 


10-5 





2-^ 


0.397 


0.731 


0.930 


0.991 


0.999 


0.9999 


1 


2-3 


0.501 


0.675 


0.921 


0.990 


0.999 


0.9999 


2 


2-4 


0.517 


0.647 


0.905 


0.988 


0.999 


0.9999 


3 


2-5 


0.524 


0.634 


0.895 


0.987 


0.999 


0.9999 



Table 4.3. Multigrid ^-cycle convergence rate (2 Gauss-Seidel Smoothers) v.s. 
jump and level 







£ 


Levels 


h 


1 


10-^ 


10"^ 


10-3 


10-4 


10-5 







0.152 


0.575 


0.904 


0.988 


0.999 


0.9999 


1 


2-3 


0.254 


0.485 


0.870 


0.984 


0.998 


0.9998 


2 


2-4 


0.269 


0.429 


0.846 


0.981 


0.998 


0.9998 


3 


2-5 


0.286 


0.403 


0.832 


0.979 


0.998 


0.9998 



Table 4.4. Multigrid y-cycle convergence rate (5 Gauss-Seidel Smoothers) v.s. 
jump and level 
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